Respiratory entrainment of the locus coeruleus modulates arousal level to avoid physical risks from external vibration

Slow rocking chairs can easily put people to sleep, while violent shaking, such as during earthquakes, may lead to rapid awakening. However, the influence of external body vibrations on arousal remains unclear. Herein, a computational model of a locus coeruleus (LC)-norepinephrine (NE) system and cardio-respiratory system were used to show that respiratory entrainment of the LC modulates arousal levels, which is an adaptation to avoid physical risks from external vibration. External vibrations of sinusoidal waves with different frequencies ranging from 0.1 to 20 [Hz] were applied to the LC based on the results of previous studies. We found that respiratory entrainment of the LC decreased the breathing rate (BR) and heart rate (HR) to maintain the HR within its normal range. Furthermore, 1:1 phase locking enhanced arousal level while phase-amplitude coupling decreased it for larger vibration stimuli. These findings suggest that respiratory entrainment of the LC might automatically modulate cardio-respiratory system homeostasis and arousal levels for performance readiness (fight/flight or freeze) to avoid physical risks from larger external vibrations.

. An illustration of the computational model for predicting arousal levels following external body vibration. Vibrational input, excitatory or inhibitory synaptic flows, CO 2 signal, mechanical signal and model outputs are depicted in blue, black, light blue, green and red, respectively. External body vibration inputs were provided by sinusoidal waves with different amplitudes and frequencies. Amygdala, LC, NTS, respiratory and cardiac centres were modelled using the Hodgkin-Huxley (HH)-type neurons, while the cardiac system was modelled using the FitzHugh-Nagumo (FHN)-type neurons. The models included different neural oscillators, wherein the default oscillatory frequencies for breathing rate (BR), heart rate (HR), and neuron groups in the brain stem were 0. 25 www.nature.com/scientificreports/ (PAC), described later in the case, without any stimulation, as reported in previous studies 20, 21 . Figure 2 shows three examples of LC-respiratory entrainment observed during the period ranging from 60 to 80 s in short-term body vibrations (condition 1). Figure 2 shows the NE-modulated conductance (representing LC activity and arousal level), the inspiration neuronal activity (representing diaphragm activity) and the expiration neuronal activity (representing abdominal muscle activity) as well as the BR and HR. , we found PAC, which modulated the phase and amplitude of LC activity by respiratory activity at 0.53 [Hz], similarly to what is observed in previous studies 28,29 . Numerical simulation of arousal level, BR, and HR elicited by short-term body vibrations. Figure 3a shows a schematic diagram of simulations performed in this study, and Fig. 3b-d show the logarithmic representation of NE-modulated conductance, BR and HR, respectively, for different frequencies with different amplitudes of vibrational input in condition 1. We noted a number of trends from the simulation results. In general, an increase in the amplitudes of the inputted sinusoidal waves tended to increase the NE-modulated conductance (arousal level), HR and BR. In a frequency range of 0.25 to 0.3 [Hz], the BR was equal to the frequency of the inputted vibration and fell within the normal range of resting BR (12-20 [bpm]) in the absolute values of the amplitude less than 12 [pA] with 1:1 phase locking, while the arousal levels were highest among all frequencies of vibration, and the HR was less than 120 [bpm] (see blue rectangles in Fig. 3b-d). At a frequency range of 0.4 to 1.3 [Hz] and an absolute value of the amplitude greater than 9 [pA], the BR was also equal to the frequency of the inputted vibration regardless of the amplitude with 1:1 phase locking. Even the BR was more than 70 [bpm], the HR almost fell within normal range of resting HR (60-100 [bpm]), while the arousal level remained high (see green rectangles in Fig. 3b-d). In a frequency range of 0.5 to 20 [Hz] and an absolute value of the amplitude lower than 6 [pA], the BR and HR almost fell within their normal ranges in the resting state, and arousal level was low with PAC (see pink rectangles in Fig. 3b-d). In a frequency range of 5-20 [Hz] and an   Fig. 3d indicate that HR significantly decreased with an increase in vibration frequency in cases with frequencies between 0.3 and 0.5 [Hz]. NE-modulated conductance could influence the neuronal activity of the amygdala via the amygdala-LC interaction modelled for reproducing fear-and stress induced activity in the amygdala. In an amplitude of − 12 [pA], the neuronal activity of the amygdala peaked at 0.2 [Hz] and decreased with higher vibrational frequencies ( Supplementary Fig. S1). Figure 4 shows NE-modulated conductance time history for different frequencies with an amplitude of − 3 [pA] in long-term body vibrations (condition 2). The arousal level changed over time, decreasing at 15 and 45 [min] while increasing at 1, 30, and 60 [min] after stimulation onset at a frequency of 0.25 [Hz]. This kind of waves are considered to be beat frequency oscillation, which is commonly observed in acoustic and optical field, where the beating represents interference between two fields of slightly different frequencies ( ω 1 and ω 2 ) and oscillates at a different frequency of ω 3 = |ω 1 − ω 2 | 30 . In our study, the red coloured line in Figure 4 corresponds to a beat frequency oscillation waveform and ω 3 was 0.00037 [Hz]  , but did not change over time (Fig. 5b). In cases of larger absolute values of vibration amplitudes more than 9 [pA], the arousal level was high in frequency ranging from 0.15 to 1.0 [Hz] and did not change over time under any frequency. In frequencies between 10 and 20 [Hz], the arousal level was lowest among all frequencies regardless of the vibration amplitudes. Further, 1:1 phase locking with PAC was observed over time, regardless of the vibration amplitudes. These simulation results in condition 2 are shown in Fig. 5 and Supplementary Fig. S2.

Discussion
This study investigated the relationship between external body vibration in daily life and arousal level. Vibration can evoke multisensory signals processed by somatosensory, vestibular and visual systems.    (Fig. 3b). In addition, neuronal activity of the amygdala peaked at a vibrational frequency of 0.2 [Hz] ( Supplementary Fig. S1). The amygdala activity influences motion sickness because the lesion of the amygdala suppress the symptoms of motion sickness 36 . These results reproduced the results of above-mentioned previous studies on the relationship between external body vibrations with particular frequencies and arousal levels induced by the vibrations. However, conflicting results have also been described, as body vibrations with frequencies of 10-20 [Hz] were reported to disturb sleepiness 8 (Fig. 4). These results demonstrate that the arousal levels change over time, which explains why a person might feel motion sickness at one time and sleepy at another at 0.25 [Hz]. This may also be one of the reasons for the contradiction that arousal level increases or decreases under external body vibration with a frequency of 0.25 [Hz] in motion sickness 6 , while a body vibration of 0.25 [Hz] induces sleepiness 9,37 . In contrast, the arousal level was low and changed little over time in cases with a frequency range of 10-20 [Hz] ( Figure 4). Our model did not determine why body vibrations with frequencies of 10-20 [Hz] disturbed sleepiness 8 . Since arousal level could be influenced by body vibration amplitudes and the entrainment of respiratory activity, further experimental studies are needed to elucidate the relationship between external body vibrations and arousal level. We investigated how arousal level changes along with external body vibrations using a computational model simulating the interaction of an LC-NE system with a cardio-respiratory system, which were modelled using the HH-or FHN-type neuron models. We noted many instances of LC-respiratory entrainment in cases with external body vibrations and found two types of pattern in LC-respiratory entrainment from the perspective of phase locking of nonlinear oscillators 17,23,29 . First, we found LC-respiratory entrainment, wherein external body vibration induced some n:m phase locking patterns with m inspiratory neuronal firings for n external vibrations, depending on the vibration amplitude. Entrainment is observed in a network including oscillator systems coupled to a pacemaker 17 , and the n:m phase locking is observed for external periodic forcing, depending on the forcing amplitude 38 . Vibration inputs with an increase in amplitude or frequency greater than 0.25 [Hz] increased LC-respiratory entrainment. The stimulus of increasing intensity or frequency can reinforce the phase response synchronisation in the neuronal population of weakly coupled neural oscillators 17,39 . In cases with frequencies between 0.4 and1. 3 [Hz] and larger absolute values of amplitudes greater than 9 [pA], 1:1 phase locking was observed (Figs. 2b and 3). We also found LC-respiratory entrainment with 1:1 phase locking in cases where vibration frequencies are close to default oscillatory ones. 1:1 phase locking was observed in cases with frequencies between 0.25 and 0.3 [Hz], which are close to the normal BR (Figs. 2a and 3). In cases with frequencies between 1.3 and 4 [Hz] and larger amplitude absolute values greater than 9 [pA], 2:1 and 4:1 phase locking were observed (Fig. 3). Second, we found LC-respiratory entrainment with PAC, in which the amplitude of the high-frequency signal is correlated with the phase of the low-frequency signal. Cross-frequency coupling of the amplitude of high frequency activity to the phase of slower oscillations in electroencephalography (EEG) recordings has been described both in animals and in humans 28,29 . In cases with a frequency range of 0.5-20 [Hz] and absolute amplitude values lower than 6 [pA] as well as in those with frequencies between 5.0 and 20 [Hz] and absolute amplitudes values greater than 9 [pA], LC-respiratory entrainment with PAC was observed (Figs. 2c and 3). The noise-induced phenomenon, stochastic resonance (SR), occurs regardless of whether the input signal is periodic or aperiodic, weak (subthreshold) or strong (suprathreshold), and regardless of whether the system is in an excitable or oscillatory regime 40 . When the subthreshold periodic signal and noise are simultaneously imposed on the oscillators, SR could manifest itself in the form of noise-enhanced phase locking 38 . The SR profile increases according to the periodic signal intensity, changing from a bell-shaped curve to a plateau, which results in the manifestation of SR without tuning of control parameters, such as the period or the amplitude of the periodic signal and the noise strength 40 . Since the model used in this study included multiple HH-type spiking neuron groups with Gaussian noise, the noise-induced SR has a positive effect on the LC-respiratory entrainment with PAC, especially in cases with frequencies between 0.5 and 20 [Hz] and smaller absolute amplitude values of less than 6 [pA].
Phase locking of oscillators is a typical phenomenon of self-organization in physical, chemical and biological systems. External periodic forcing above a threshold induces an n:m phase locking pattern with m firings for n forcing periods, depending on the forcing amplitude 38 . As shown in Fig. 3b-d, we noted different types of LCrespiratory entrainment with n:m phase locking. With an increase in the frequencies of external body vibrations from 0.4 to 20 [Hz], the mode of LC-respiratory entrainment shifted from 1:1 to 2:1 and 4:1 phase locking and then to PAC, with the BR significantly decreasing in case when the absolute amplitude values were greater than 9 [pA], which is shown as orange lines in Fig. 3(c). Fig. 6 shows an example of mode shifts from LC-respiratory entrainment with 1:1 phase locking to PAC in cases with an amplitude of − 12 [pA] and frequencies ranging from 1.1 to 10 [Hz] in condition 1, and indicates that phase locking between LC activity and inspiratory activity was found in four cases, and BR decreased with an increase in vibration frequency. Fig. 7 shows comparisons of neuronal activity of LC-respiratory, retrotrapezoid nucleus (RTN), inhibitory IE, excitatory IE, and NTS (Baro) in mode shifts from LC-respiratory entrainment with 1:1 phase locking to PAC in the same cases shown in Fig. 6. At a frequency of 1.1 [Hz], activities of RTN, NTS (Baro), excitatory IE and inhibitory IE were entrained by both LC activity and respiratory activity with 1:1 phase locking. However, at frequencies of 1.7 and 3 [Hz], activities of RTN and inhibitory IE were entrained by LC activity, while those of NTS (Baro) and excitatory IE were entrained by respiratory activity. At a frequency of 10 [Hz], activities of RTN, NTS (Baro), excitatory IE and inhibitory IE were entrained by respiratory activity, and PAC was generated. The RTN directly received excitatory synaptic signals from the LC, while the NTS (Baro) received mechanical feedback by pulmonary stretch receptors located in the lungs and transmitted information on the lung volume, that is, the BR (see Fig. 8). Therefore, the RTN was strongly influenced by LC activity, while the NTS (Baro) was strongly influenced by the respiratory activity. The excitatory IE received excitatory signals from ramp-I, while the inhibitory IE received excitatory signals from both ramp-I and late-E, which in turn directly received excitatory signals from the RTN. Thus, the inhibitory IE was influenced by RTN as compared to the excitatory IE, especially in cases with active expiration (shown as abdominal muscle activity in red lines, the 1st row of Fig. 7) to which quantal acceleration of late-E activity leads 41 .
Respiratory entrainment has a significant role on the mode shift. Respiration entrained the activity of each neuronal group, such as RTN, NTS (Baro), excitatory IE and inhibitory IE, in addition to modulating the default oscillatory frequency of each group to BR. In a case with a frequency of 1. . The reason of the decrease in HR is due to the differences in neuronal activity between RVLM and AMB (see Supplementary Fig. S3). In the cases with a frequency range of 0.3-0.5 [Hz], 1:1 phase locking enhanced the neuronal activity of ramp-I and late-E in the respiratory centre with an increase in the BR and modulated the default oscillatory frequencies of RTN, NTS (Baro), excitatory IE, Inhibitory IE, RVLM and AMB to those of the external vibrations, where active expiration induced an increase in the activity of late-E and then increased the activity of AMB to a greater extent than that of RVLM. As a result, the increase in vibration frequency from 0.3 to 0. 5  , even at larger frequencies. As mentioned previously, the noise-induced SR has a positive effect on the LC-respiratory entrainment with PAC and may stabilise the BR and HR within their normal ranges in the resting state. Previous studies have demonstrated the potential benefit of using the noise-induced SR, in which an auditory noise significantly improves the balance control of participants with lower balance control ability 46 . Therefore, LC-respiratory entrainment might have a role in maintaining homeostasis of the cardio-respiratory system.   48 . An increase in the resting BR and HR serves as alternative indices for measuring lifethreatening risks, because hypovolemia induced by bleeding causes an increase in sympathetic nerve activity, leading to an increase in the BR and HR. In addition, BR increased up to 60 [bpm] during exercise using muscle activity such as running and high intensity internal training 49 . Since the vibration stimulation can increase muscle activity to stabilize the human posture, it may increase BR to higher values.  Previous studies have demonstrated that the mode shift from 1:1 to 2:1 phase locking has some benefits for enhancing the auditory sensitivity and selectivity of male mosquitoes who detect flying females 50 as well as for enhancing the salience of weak signals for high-power lasers in a technique known as injectionlocking 18,50 . Therefore, our finding on the mode shift in cardio-respiratory system modulation might represent a biological benefit of entrainment, although further neurophysiological studies are needed to confirm it.
Although previous research has shown that the external body vibrations within a frequency range of 0.1-20 [Hz] investigated in our study may pose physical risks, we did not find any study that directly investigated the effects of the earthquakes on human body. The ground shaking due to earthquakes causes buildings to vibrate at frequencies in a range of approximately 0.1 to 30 [Hz] 51 , and it is estimated that the human body experiences  41 . The closed-loop respiratory system includes two major feedback pathways from the lungs to the respiratory CPG: mechanical feedback by pulmonary stretch receptors (PSR) and chemical feedback by peripheral chemoreceptors as well as central chemoreceptors from retrotrapezoid nucleus (RTN) neurons. The cardiac centre comprises pathways that modulate heart activity: sympathetic outflow with norepinephrine (NE) modulation from the rostral ventrolateral medulla (RVLM) and parasympathetic outflow via acetylcholine (Ach) modulation from the ambiguous nucleus (AMB). Inspiratoryto-expiratory phase-spanning neurons (IE) in the pons combines the activity of the respiratory and cardiac centres. Neurons were modelled using Hodgkin-Huxley (HH)-type spiking neurons. Vibrational input is directly provided to LC. Heart rate (HR), breathing rate (BR) and NE-modulated conductance representing arousal level are outputted from cardiac system, lung and LC, respectively. cVRG: caudal ventral respiratory group.  15,16 , while the vibrational sinusoidal stimulations to vestibular sensory organs for inducing responses of caudal medullary raphe neurons were performed for forward and backward rotations (pitch) around the lateral direction 54 . However, our model cannot set the direction of body vibration to the longitudinal, lateral and vertical directions because it did not include caudal medullary raphe neurons or other related neurons. Further experimental and computational studies are necessary to investigate effects of vibration direction on the arousal level. Finally, the integrated model of the LC-NE system and cardio-respiratory system included two neuronal oscillator groups, except the respiratory system, that is, the cardiac system and neuron groups in the brainstem, which have default oscillatory frequency ranges of 1-1.6 [Hz], and 4-10 [Hz], respectively. The system could have nonlinear oscillator ensemble to create complex phase locking patterns 24 . However, we could not determine how default oscillatory frequency ranges of the two neuronal oscillators are related to the LC-respiratory entrainment we found in this study. This might be because the respiratory system can have a much stronger entrainment effect than the other neuronal oscillators. Future studies should investigate the mechanism of LC-respiratory entrainment using a computational model including more neuronal oscillators simulating various ones found throughout the brain and body.

Methods
Coupling of the LC-NE and cardio-respiratory systems. Figure 8 illustrates the model used in this study. We modelled the amygdala, RTN, RVLM, caudal ventrolateral medulla (CVLM) and AMB using 20, 40, 30, 30 and 30 HH-type spiking neurons, respectively. Baroreceptor-related neurons in the NTS (NTS Baro) and IE include 20 and 30 sets of excitatory and inhibitory HH-type neurons, respectively. Chemoreceptors in the NTS (NTS Chemo) include two groups of 30 HH-type neurons: one receives CO 2 signals from the lungs and excites or inhibits the amygdala 55,56 , and the other comprises second-order peripheral chemoreceptors, which have a direct projection to the RTN, preI/I (pre-BötC) and RVLM 57 . Since the number of HH-type spiking neurons may influence the simulation results 58 , we performed parametric simulations using different number of HH-type spiking neurons. The default numbers of HH-type spiking neurons for the amygdala, RTN, RVLM, CVLM, AMB, NTS Baro, NTS chemo and IE were set to 20, 40, 30, 30, 30, 20, 30 and 30, respectively. We simulated 0.5, 1.5 and 2 times of the default numbers of HH-type spiking neurons, for instance, in case of 2 times, they were set to 40, 80, 60, 60, 60, 40, 60 and 60, respectively. The parametric simulation results demonstrate that smaller numbers of HH-neurons had different values in BR and HR in comparison with the default numbers ( Supplementary Fig. S4).
Each neuron can be fundamentally represented using the following equations.
where C is the membrane capacitance as C = 36 [pF], while g K , g Na and g L are the peak conductances of potassium and sodium as well as the leak conductance (Leak), respectively. g K = 250 [nS], g Na = 400 [nS] and g L = 6 [nS], while E K , E Na and E L represent the reversal potentials of potassium, sodium and Leak, respectively. The constants were determined based on 59 , as follows: www.nature.com/scientificreports/ The time constants were set as follows: where τ mK and τ hNa are in [ms]. m Na was assumed to be m ∞N a , and the time constant was assumed to be τ mNa = 0 based on a previous result 60 . ξ(t) is Gaussian noise, and D is noise strength, with D = 20 [pA]. Equation (1) was updated using the Euler method. In this study, the time step was set as 1.0e−4 [s]. The LC was modelled using a 15 × 15 excitatory neural network with an HH-type equation, as follows: where I E syn is the excitatory synaptic input, and I I NE represents the NE interaction, as described later. I gap is the electrical current based on the gap-junction from the eight neighbouring neurons and is computed as . ϕ is the modulating term reflecting CO 2 concentration (pCO 2 ) in the blood, and was described according to a previous study 61 61 used an HH-type model that included Ca 2+ channels, we used one without them for simplicity. NE emitted by the excitatory firing of LC neurons is diffused to cognitive neurons and modulates the activation gain g NE . The modulatory dynamics for g NE are described using the second-order delay system, as follows: where t j is the jth firing spike, and g NE corresponds to the NE-modulated conductance released into the surroundings.
. The rise and decay time constants of τ 0 and τ 1 are 100 [ms] and 300 [ms], respectively, as previously reported 62 .
The LC in the pons supplies NE input to many regions of the cortex and the respiratory centre. The diffusion of NE in the cortex is associated with arousal; an increase in NE is accompanied by an increase of arousal level. Chemosensitivity in LC neurons is activated by an increase in CO 2 concentration in the blood, as the LC receives excitatory information from the respiratory CPG (central pattern generator). We found that the controller gain of the chemical feedback from the lung to the LC can influence the time-to-respond for an emotional fear-or surprise-related stimulus 63 . LC neurons also receive excitatory input from preI/I (pre-BötC) neurons, as previously shown in adult mice 11 . Since LC neuron activity is synchronized with inspiratory bursts, we modelled the connectivity from preI/I to LC with all-to-all connections. However, there is no evidence of their synaptic strengths. We set their synaptic strengths to 0.1 based on synaptic strengths for the other connectivity in respiratory CPG. In addition, the direct pathway from the LC to the RTN/pFRG (parafacial respiratory group) induces active inspiration and expiration 64 . Therefore, NE diffusion, that is, the activity of LC, indirectly influences the control of breathing. We included the amygdala-LC interaction in the model because fear-and stress-induced activity in the amygdala influences sensory brain regions through connections with the LC 65 . Furthermore, increased tonic activity in LC neurons induces anxiety-like and aversive behaviour 66 . However, we did not include the cortical and subcortical controls of breathing associated with emotion and cognition. Although the amygdala-LC interaction includes GABAergic inhibitory connections 67 and corticotropin-releasing factor, which has an excitatory effect on LC neurons 68 and adjusts the activity as well as reactivity of the LC-NE system 69 , we only included the excitatory interactions between the amygdala and LC in our model.
The respiratory centre and respiratory system were developed based on a mathematical model of the closedloop system for the control of breathing proposed by 41 . Breathing or lung ventilation represents the exchange of air between the lungs and the surrounding environment, which takes place according to the rhythmic contraction  www.nature.com/scientificreports/ of the inspiration and expiration neurons. The firing activities of the phrenic and abdominal motor neurons control the inspiration and expiration neurons, respectively, representing major outputs of the brainstem respiratory CPG, which generates respiratory oscillations. The closed-loop respiratory system includes two major feedback pathways from the lungs to the respiratory CPG, namely mechanical feedback and chemical feedback. Mechanical feedback is provided by pulmonary stretch receptors (PSR) located in the lungs, which transmit information regarding lung volume to the brainstem via the vagus nerve. Chemical feedback is assumed to be provided by changes in the peripheral CO2 concentration as a signal from the peripheral chemoreceptors to NTS Chemo 57 , which is sensitive to the levels of oxygen (O 2 ) and CO 2 in the blood and promotes RTN activity, thereby causing an increase in HR and BR (see Fig. 8). We adjusted the NTS Chemo parameters so that active expiration was not induced by the activity of abdominal muscles during normal breathing; however, active expiration was induced during hypercapnia according to previous results 57 . The arterial CO2 partial pressure of around 40 mmHg was provided to central chemoreceptors of the RTN during normal breathing. The central chemoreceptors of the RTN, which are highly sensitive to the CO 2 concentration in the brainstem, and thus modulate respiratory CPG in a CO 2 -dependent fashion. The parameters of the respiratory centre and system were originally adjusted to the adult rat respiratory system 41 . We replaced mechanical system of the lung with a mathematical model of respiratory muscles, chest wall and lungs based on adult human subject data previously proposed by 70 . Then, we modified the parameters to adjust the adult human respiratory system as listed in Table 1.
The cardiac centre was developed based on the pathways via which the central nucleus of the amygdala (CeA) influences blood pressure during mental stress or anxiety, as proposed by 72 , and autonomic chronotropic control of the heart, as previously depicted 73 . During stressful conditions, the CeA may inhibit baroreceptive neurons in the NTS. This inhibition might switch off inhibitory inputs from the CVLM to the RVLM, which could activate neurons in the RVLM, leading to an increase in sympathetic outflow with NE modulation and an increase in blood pressure and HR. The AMB also receives excitatory information from the NTS, which increases parasympathetic outflow via acetylcholine (Ach) modulation of heart activity, thus decreasing HR. The cardiac system was represented using a cardiac muscle activation model with electrical activity based on the FitzHugh-Nagumo (FHN) neuron model 74 , as well as the neuromodulation of NE and Ach, which are neurotransmitters in sympathetic and parasympathetic nerves, respectively. The electrophysiological model of cardiac muscle activity is described by the following diffusion equation: where v is the membrane potential of cardiac muscle cells, and w is the recovery variable. c 1S = G c c 1S , c 2S = G c c 2S , c 1S = 1 , and c 2S = 0.22 . α , b and d are parameters that define the shape of v. We hypothesised that G c and b would be determined by the neuromodulation of NE and Ach, as shown in the following Eqs. (20,21), based on Hill's Eq. 75 .
The HR increases with the value of parameter b and decreases with the value of b. H(q Ach ) significantly contributes to the adjustment of HR. The average HR ranges from 64 to 69 [bpm] in adult men aged 20-60 years 26 , H(q Ach ) = 1 − 1/ 1 + q Ach mid /q Ach α Ach .  27 . We set p NE mid = 350.0 , q Ach mid = 20.0 , α NE = 10.0 , and α Ach = 5.0 to reproduce the average values of BR and HR.
In this study, we introduced inspiratory-to-expiratory phase-spanning neurons (IE), including excitatory and inhibitory neurons, into the pons to combine the activity of the respiratory centre with that of the cardiac centre. The excitatory IE was proposed by 57,76 to receive an excitatory signal from ramp-I of the respiratory centre and project to the RVLM, providing a pontine-dependent inspiratory modulation of thoracic sympathetic nerves according to a suggestion of respiratory-sympathetic coupling based on previous experimental observation 77 . Further, the role of the excitatory IE was confirmed with additional computational and experimental studies 78 . The inhibitory IE was described to receive excitatory signals from the respiratory centre, projecting into the AMB and providing respiratory-related modulation of parasympathetic nerves based on recent neurophysiological studies 79,80 . Here, we added an excitatory projection from ramp-I neurons in the respiratory centre to the inhibitory IE in order to control how the flow of information from the NTS to the AMB is cut off during inspiration, as per 73 , as well as an excitatory projection from late-E neurons to inhibitory IE to represent hypoxia-induced transient tachycardia, as reported by 81 .
This model has several parameters; we mostly determined the parameters based on previous studies 41, [60][61][62]70,71 , However, we used four parameters p NE mid , q Ach mid , α NE and α Ach to adjust the BR and HR to be 16 and 68 [bpm], respectively, which are almost the same as those in healthy adult males in the resting state, by approximately five times of trial and error.
Experimental design. We used the developed computational model for autonomic breathing simulation under the following two types of external body vibration conditions (Fig. 1) , hereafter referred to as extraction points. As mentioned previously, vibrational inputs to mammalian bodies are known to stimulate vestibular sensory organs and then induce LC activity that is physically synchronised with the original vibrational inputs 15 . However, how the amplitudes of vibrational inputs are transmitted to the LC is unknown. Therefore, we used current values inputted to HH-type neurons to represent the amplitudes of the stimulations inputted to the LC. All amplitudes represent excitatory stimulation, which increased with the absolute value.
We calculated the BR, HR and NE-modulated conductance from the LC to the cortex, which represent the arousal level, activity of the inspiration and expiration neurons, as well as neuronal activity of amygdala, NTS (Baro), RTN, IE, RVLM and AMB as model outputs. Activity of the inspiration and expiration neurons ranged from 0 to 1, while neuronal activity of amygdala, NTS (Baro), RTN, IE, RVLM and AMB was calculated as neuronal spike counts. The BR, HR and NE-modulated conductance were calculated as averaged values during the period ranging from 60 to 80 [s] in the autonomic breathing simulations for condition 1 and during the period of 20 [s] before each extraction point for condition 2. We investigated the relationship between external body vibration input and model outputs of arousal level, BR, and HR using the obtained simulation results. The model was implemented using C++ language, and all simulations were performed using a computer with an Intel Xeon Gold 6242 (16C/2.8Ghz) and 384 GB Memory.

Data availibility
The source codes used to generate the results in this paper are available at https://osf.io/bzdtf/?view. www.nature.com/scientificreports/